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In previous articles [J. Chem. Phys. 121 4501 (2004), J. Chem. Phys. 124 034115 (2006), J. 
Chem. Phys. 124 034116 (2006)] a bipolar counter-propagating wave decomposition, <]/ = >]/+ + <]/-, 
was presented for stationary states ^ of the one-dimensional Schrodinger equation, such that the 
components $± approach their semiclassical WKB analogs in the large action limit. The corre- 
sponding bipolar quantum trajectories are classical-like and well-behaved, even when has many 
nodes, or is wildly oscillatory. In this paper, the method is generalized for multisurface scatter- 
ing applications, and applied to several benchmark problems. A natural connection is established 
between intersurface transitions and (+<-►—) transitions. 
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I. INTRODUCTION 



This paper is the fourth in a series^™ investi- 
gating the use of "counter-propagating wave meth- 
ods" (CPWMs)ii2i2i^ for solving the time-independent 
Schrodinger equation exactly. The basic idea is to decom- 
pose the stationary wavefunction, 5", into a two-term, or 
"bipolar" form, 



(1) 



such that the \&+ and "f- represent wave "components" 
moving in opposite directions. Although Eq. ([TJ) is exact, 
the particular bipolar decomposition used is chosen to 
correspond to analogous approximate semiclassical wave 
components'^ 9 - For smooth potentials, the latter are 
known to exhibit smooth and slowly-varying field func- 
tions throughout the interaction region. The exact *5f± 
of Eq. |T]) must behave similarly, at least in the classical 
limit, and in any case lead to field functions that are much 
better behaved than for W itself. For stationary scatter- 
ing states, for instance, ^ necessarily exhibits oscillatory 
interference in one or more asymptotic regions where the 
interaction potential V is flat, whereas the asymptotic 
$4. and \&_ behave as plane waves. Thus, interference 
manifests not in the individual components, but arises 
naturally from their linear superposition. 

Apart from certain conceptual and theoret- 
ical advantages, the above picture is particu- 
larly relevant for quantum trajectory methods 
(QTMs) ,milil2ii2ililLlilLim2iU-i.e. , trajectory- 
based numerical techniques for performing exact 
quantum dynamics calculations, in a manner similar to 
classical simulations; 21 QTMs that are based on stan- 
dard Bohmian mechanic o 22 ' 23 ' 24 ' 25 ' 26 ' 27 use a single-term 
or "unipolar" representation of the wavefunction, from 
which the quantum trajectory evolution is determined. 
However, the Bohmian time evolution equations are 
nonlinear, which can lead to radically different QTM 



behavior when applied to the individual bipolar compo- 
nents of Eq. |T]), rather than to W itself. In particular, 
the smooth field functions of the ^f± obviate the "node 
problem" — i.e. the near-singularities in the quantum 
potential, Q, leading to numerical instabilities in the 
quantum trajectory evolution when there is substantial 
interference! 5 ' 6 ' 15 ' 28 ' 29 ' 30 A more detailed discussion may 
be found in the previous articles cited above. 

For one-dimensional (ID) stationary scattering state 
calculations, the bipolar CPWM trajectories are actu- 
ally classical, in a certain generalized sense (Sec. |TT] 
and Ref. |J). Quantum effects arise not from a quan- 
tum potential, Q, impacting trajectory evolution, but 
rather through dynamical coupling between the two 
components, and \f r _. This coupling, which is in 
essence proportional to the interaction potential, induces 
(+ <-> — ) transitions — i.e., scattering from the , J+ in- 
cident/transmitted wave to the reflected wave, and 
vice- versa. In many respects, the situation is reminiscent 
of traditional multisurface scattering theory, in which an 
off-diagonal diabatic coupling potential induces a dynam- 
ical transition from one diabatic state to another. In- 
deed, one somewhat compelling conclusion of the present 
work — which deals with a bipolar CPWM treatment of 
multisurface dynamics — is that in some sense both types 
of transitions may be regarded as different aspects of the 
same underlying phenomenon. This "unification" may 
result in an interesting cross-fertilization of ideas, e.g. in 
the area of classical trajectory surface hopping (TSH) 
methodologies, originally designed by Tully i 31 ' 32 

We are by no means the first to apply QTMs to the 
dynamics of electronic nonadiabatic collisions; the first 
papers to do so were written several years ago by Wyatt 
and coworkers i 14 ' 17 It should be noted that Wyatt's treat- 
ment is exact — at least formally — although approximate 
and/or classical versions have also been developed ! 33 ' 34 ' 35 
These methods are often compared to TSH, with which 
they (and the present approach) have important differ- 
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ences. In TSH, a swarm of independent classical trajec- 
tories evolve along a given potential energy surface; for 
each trajectory, a decision is made whether to "hop" to a 
different surface, based upon a transition probability, ob- 
tained by integrating coupled equations for the electronic 
amplitudes. In contrast, the above multisurface meth- 
ods do not involve any trajectory hopping (the number 
of trajectories evolving on each electronic surface is, in 
fact, conserved) but transfer density and phase informa- 
tion from one state to another in a continuous manner. 
Note that — like its single surface counterpart — the mul- 
tisurface QTM of Wyatt et al. is based on a unipolar 
representation of the wavefunction, and is therefore also 
subject to numerical instabilities due to interferences. 
Moreover, in the multisurface case, interference arises not 
only from waves propagating in opposite directions on the 
same surface, but also from waves transferring flux from 
one surface to another. 

It should also be stated that in the multisurface con- 
text, the idea of applying a bipolar decomposition has 
been previously considered. In particular, M. Alexan- 
der and coworkers adopted such a scheme in the exact 
quantum solution of the close-coupling equations using 
log-derivative propagation ) 36 ' 37 although their choice of 
Eq. (TT]) decomposition does not avoid oscillatory field 
functions, and is therefore not so useful for QTMs. D. V. 
Shalashilin also used a bipolar decomposition to solve the 
close-coupling equations (for inelastic scattering applica- 
tions), albeit only as a semiclassical approximation ] 38 ' 39 

The remainder of this paper is organized as fol- 
lows. Sec. |n] discusses theoretical and algorithmic de- 
velopments, both for single surface dynamics calcula- 
tions (Sec. Ill Ap and the multisurface case (Sees. Ill Bl 
and III Cp . Note that the bipolar CPWM algorithms used 
here require neither complex scalin g 40 ' 41 ' 42 nor absorb- 
ing potential o 43 ' 44 ' 45 ' 46 — a decided advantage over other 
quantum scattering methods. Moreover, the scaling of 
computational (CPU) effort is linear with the grid size, 
N, rather than proportional to N 3 . Results are presented 
in Sec. IIIII For single surface applications (Sec. IIII Ap . 
in addition to a test suite of wide-ranging applications 
considered previously^ we investigate the ability of the 
method to compute scattering quantities to extremely 
high relative accuracy, or precision. Several benchmark 
two-surface applications are considered in Sec. IIII Bl in- 
cluding two of the Tully models. 3 - Concluding remarks 
are given in Sec. [IV] 



II. THEORY AND ALGORITHM 
DEVELOPMENT 

A. Single Surface Dynamics 

In several previous articlesj 2 ^ 3 -'^ an extremely accurate, 
efficient, and robust bipolar CPWM numerical algorithm 
was developed for computing stationary scattering states 
of ID single-surface systems. A brief summary is pre- 



sented here; further details can be found in the above- 
referenced articles. 

The algorithm is a time-dependent relax- 
ation method, for which the initial wavefunc- 
tion is a left-incident/transmitted wave only, i.e. 
^(x,t = 0) = %+(x,t = 0) and *_(x,t = 0) = 0. 
Over time, a reflected wave ty-(x,t) comes into being 
through interaction region coupling due to the potential 
energy, and eventually ^(x,t) — ^ + (x,t) + \P _ (a;, £) 
relaxes to the true stationary eigenstate of desired 
energy, E, and left-incident boundary conditions. The 
solution components ^±(2;) behave as (right/left) 
traveling plane waves in both asymptotes, except that 
lirn^oo ^_(x) = 0. The asymptotic square amplitudes, 
lim^oo |^I' + (j:)| 2 and lim^-oo |^_(a;)| 2 , are directly 
related to transmission and reflection probabilities, 
respectively [Eq. ((H) ] . Through the interaction region, 
the solution ^±(a;)'s are exact quantum analogues of a 
type of semiclassical WKB approximation resulting from 
the "generalized Froman" (F) approach^! 8 - Amplitude 
and phase functions for the components are expected 
to be smooth and slowly varying, unlike those for the 
solution itself. 

The "generalized" aspect of the above methodology 
implies that we are formally allowed to specify the "clas- 
sical" trajectories as we wish, using an effective poten- 
tial V cii (x) of our own choosing, rather than the ac- 
tual potential V{x). In practice, the method requires 
that lim^ioo [V(x) — V eS (x)] = 0, and works best 
when V cS (x) is smooth and monotonic between the two 
asymptotes. This ensures that all trajectories are de- 
void of turning points, and that there is no asymptotic 
coupling between the two ^± components. Note that 
the 1 &+(x) trajectories move to the right with velocity 
v(x) = y/2 [E — V cS (x)] /m, and the trajectories 
move to the left with equal and opposite velocity — v(x). 
As per other QTM's, the trajectories give rise to discrete 
moving grids, which carry local phase and amplitude in- 
formation for the appropriate \l/± component. There is 
no quantum potential 6 -^ contributing to the dynamics; 
all trajectories are (generalized) classical, and can be de- 
termined a priori. Instead, quantum effects manifest via 
ty± coupling in the interaction region — as evidenced by 
the time evolution equations (Ref. |j Eq. (16), or Eq. ([7]) 
of this article with / = 1). The coupling is essentially 
proportional to the interaction potential, and can be said 
to induce transitions from one \l/± component to the 
other. 

In the previous papers, various numerical innovations 
were introduced to improve performance and stability of 
the trajectory-based time propagation of \l/±, designed 
to deal effectively with the unusual mixed boundary con- 
ditions [Eq. (fT5|) ]. First, since the upper (+) and lower 
(— ) grids are moving in opposite directions, they can- 
not be commensurate at all times. However, a scheme 
was introduced whereby these two grids are commensu- 
rate at all times t equal to a multiple of the "shift time," 
^shift — the time it takes one grid point to travel to the lo- 
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cation vacated by its nearest neighbor. A constraint on 
ishift is that it be a multiple of the propagation time step 
size, A. Second, since grids are generally incommensu- 
rate, interpolation must be used to evaluate the coupling 
contribution — a procedure found to yield far better re- 
sults when applied to a polar (amplitude/phase) decom- 
position of the \&±'s, rather than directly to the ty±'s 
themselves. A "plane wave propagator" (PWP) approach 
was also employed, which effectively treats the uncoupled 
part of the time evolution equations exactly, and the cou- 
pled part using first-order Euler time integration. 



Unfortunately, the PWP idea is found to be generally 
incompatible with more efficient time integrators, such as 
fourth-order Runge-Kutta (adaptive and non- adaptive), 
used here. Note that for symmetric potential systems, 
it is possible to recast the time evolution equations in 
such a way that the PWP contribution vanishes (by in- 
troducing a time-evolving phase into the definition of 
the wavefunction) . In order to ascertain the extent to 
which PWP actually improves performance, some pre- 
liminary numerical tests were conducted for the Eckart 
A system (Sec. with E = Vq. In particular, an effi- 
ciency comparison was made between the original PWP 
algorithm and the "phase-modified" version — both using 
first-order Euler propagation for the non-PWP contri- 
bution. For an identical set of numerical parameters, 
both methods performed about equally well, vis-a-vis 
the accuracy of the computed reflection and transmis- 
sion probabilities — although the phase-modified version 
executed roughly 30% faster, owing to less required com- 
putation. As a second, more realistic test, the phase- 
modified algorithm with fourth-order Runge-Kutta was 
compared to a straight fourth-order Runge-Kutta treat- 
ment of all terms in the original evolution equations (in- 
cluding PWP terms). Both methods again performed 
roughly equally well with respect to computed accuracy, 
with slight performance differences depending on the de- 
sired level of accuracy. Moreover, both methods are far 
superior to the first-order Euler methods, enabling A val- 
ues that are orders of magnitude larger, without com- 
promising numerical accuracy. Only the straight Runge- 
Kutta method can be generalized for asymmetric poten- 
tial systems, however, and is introduced here as the gen- 
eral method of choice. 



We also introduce an improved method for dealing with 
extremal trajectories that stray outside the range of in- 
terpolation: instead of ignoring the coupling contribu- 
tion altogether^ we now simply use values for the polar 
field quantities obtained from the nearest extremal grid 
point — although it should be emphasized that the true 
coupling contribution is vanishingly small in the asymp- 
totic limit. Also, even within the range of interpolation, 
the interpolated density value may become slightly neg- 
ative when the density absolute values are very small; in 
that case, we now reset the interpolated value to zero. 



B. Multisurface Dynamics: Theory 

We now generalize the theory for the ID multisurface 
case. Note that only a brief summary is provided here, 
as the full derivation follows closely that of Refs. 3 and 4, 
which should be consulted for further details. 

Let / be the number of electronic states considered, 
with no restrictions on intcrsurfacc coupling. A diabatic- 
like time-independent matrix Schrodinger equation is 
presumed, of the form 



(2) 



where {^i, ^2, •••^ r /} comprise the vector components 
(associated with each of the / diabatic states) of the nu- 
clear wavefunction, 'J, and 
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are the components of the / x / Hamiltonian operator 
matrix, H, with i < / and j < / labeling diabatic states. 

The Vij(x) = Vj : i(x) are the diabatic potential energy 
curves, with the i ^ j case denoting the coupling poten- 
tials. In order to ensure that coupling vanishes in the 
asympotic limits (required to obtain asymptotic scatter- 
ing waves with correct boundary conditions) we must 
have ]im x —>± 00 Vi^ : j(x) = 0. However, the asymptotic 
values for the diagonal potentials, Va(x), are allowed 
to be completely arbitrary, and in particular, need not 
be symmetric. Left and right asympotic values are de- 
noted V iL = lim^-oo Vi,i(x) and V lR = lim^oo Vi,i(x), 
respectively. For purposes of generating trajectories 
for motion on the i'th diabatic state, a suitable family 
of / effective potentials, Vf s (x), are introduced, such 
that lim a! _,± 0O [Vi,i(x) — V^ e (&)] = (to ensure against 
asymptotic coupling), and V i cS (x) < max(ViR, Vn) for 
all x (to avoid turning points, as discussed in Sec. Ill Al) . 

Assuming that each component of is decomposed 
into its own bipolar expansion, 



%(x) = y i+ (x) + y i _(x) 1 



(4) 



we clearly need 2/ independent differential equations 
[and associated boundary conditions (Sec. Ill C|) ]. to 
uniquely specify the solution decomposition of Eq. |4|. 
One half of these are already provided in Eqs. ^ and 
([3]). For the remaining / equations, it is natural to ap- 
ply the generalized F approach (Sec. Ill A|) in component- 
wise fashion. In particular, for a single-surface system, 
the generalized F approach 8 provides the following as the 
second (after the usual Schrodinger equation) indepen- 
dent differential condition, 



v imv , 
2v n 



(5) 



(where the prime denotes spatial differentiation), as dis- 
cussed in more detail in Ref. i [Eq. (11)]. For the mul- 
tisurface generalization, we simply apply Eq. (JSJ) to each 
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diabatic state separately, to obtain 



2vi 



(6) 



where Vi(x) — y2[E — V? s (x)] /m are the trajectory 

velocities for the z'th diabatic state (as determined by the 
effective potential V° s ). Note that Eq. ([6]) does not de- 
pend at all on the Vi^j as appropriate, i.e. all intersurface 
coupling should arise through the Schrodinger Eq. ([2]), it- 
self, if (+<->—) transitions for single- and multi-surface 
applications are to be treated equally. 

By combining Eq. ^ with Eq. we can derive ex- 
pressions for the first spatial derivatives of each of the 
2/ bipolar CPWM components, i.e. the \&i±' , directly in 
terms of the undifferentiated component quantities. This 
yields somewhat complicated results, analogous to Ref. [1 
Eq. (12), which are excluded here for the sake of brevity. 
By following the procedure described in detail in Ref. 3 
and Ref. |J [basically, constructing a convective term to 
get rid of first-order spatial derivatives of the wavefunc- 
tion in the hydrodynamic frame, and introducing an ex- 
plicit time dependence via d 1 $±/dt = — (i/h)E^f±\, we 
are then led to the following coupled hydrodynamic (La- 
grangian) time evolution equations: 




[E-Vu-VfS + a 



where, 



2m 



V'- 



16 \ E — Vi 



cir 
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4 £- Vf s 



(8) 



For pedagogical purposes, we also consider the asympot- 
ically symmetric special case with Vn, = Vm = and 
Vf s (x) = for all i: 



dt 



f 



E* i± -J2ViA*. 



3 + 



*4 



(9) 



Note: in Eq. (J5|), j ranges over all values, including j = i. 

The time evolution equations of the preceding 
paragraph — best exemplified by Eq. ((9]) — offer a coher- 
ent, unified picture of scattering theory that places both 
intersurface transitions and (+ «-* — ) transitions on a 
near-equal footing. Note first that coupling of all types 
vanishes in both asymptotic limits, resulting in uncou- 
pled, scattering plane wave dynamics in these limits for 
all \&i±, as desired. In particular, i &i±(x — ► ±oo, t) is 



a plane wave propagating to the (right/left) with uni- 
form speed ViL/n — ^2 [E — V^L/ii] l m - ^ i s the var- 
ious potential interactions, operating in the interaction 
region, that induce transitions of all kinds among the 
2f bipolar CPWM components, ^fi±. The off-diagonal 
potentials Vi^j are responsible for transitions between 
diabatic states i and j, whereas the diagonal potentials 
Vi t i induce transitions from reactive (transmitted) ^,-4. 
to non-reactive (reflected) components (and vice- 
versa) . Apart from this distinction, both types of transi- 
tions manifest similarly in the equations of motion. 

As in the previous workj^ the time evolution equa- 
tions above give rise to some elegant flux properties, 
which again treat both types of transitions (or proba- 
bility flow) on a near-equal footing. Let ji± = ^iPi± 
be the flux for component i &i±, defined in terms of the 
generalized trajectory velocities Vi, and the component 
probability densities, pi± = l^iil 2 . We also find it con- 
venient to introduce new composite labels, a — (i, ±) and 
j3 = (j, ±), to label individual *ffi± components. Note 
that the ± values are independent for a and (3, so that 
each of these two component labels can take on 2/ dis- 
tinct values. Using Eq. ([7]), and transforming to Eulerian 
(partial) time derivatives, we obtain the following flux re- 
lation: 



dp a _ . 1 
dt ~ Ja 



0^a 



dt 



(10) 



where 



dp" 



-0 



dt ~ n [Vh > 
-0 , 



J i,3 



{Vf + Ci)]lm[^ a *^ ], (11) 



and dp"^ /dt represents the rate of probability density 
flow from component (3 to a. 

Equations (fTU)) and (fTTj) above emphasize the essen- 
tial similarities between the two types of transitions — 
particularly for the Vn, — Vm — case, for which the 
8ij contribution in Eq. Ijlip vanishes. More importantly, 

however, we find that dp^^/dt = —dp^ a /dt, which 
implies the nontrivial combined continuity relation.;^ 



dt 



(12) 



Equation (JT2]) further implies that the total probabil- 
ity integrated over all 2/ components is conserved over 
time — a desirable but nontrivial result, given that 
and ^>i^ are not true "components" (orthogonal comple- 
ments) in the way that \&j and <f j^i are. 

The above discussion regarding time evolution proper- 
ties pertains to all times, t — not just the t — > 00 limit, 
where VP approaches the exact stationary solution. In 
the latter limit, however, it is clear that Eq. (fT2"|) must 
be zero, so that the stationary solution must satisfy 
j+ = —j-', where j± = X)f=ii«± represents the to- 
tal flux moving in the (right/left) directions. For the 
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Vil = ViR = case in particular, this implies that 
p+' = p-' , where p± = ^2 i=1 Pi±- In other words, the 
sum of all + densities is equal to the sum of all — den- 
sities, apart from a constant — a nice generalization of 
Ref.iEq. (15). 

C. Multisurface Dynamics: Numerical Algorithm 

To implement the above time evolution equations nu- 
merically, we must first discuss boundary conditions and 
initial value conditions. Without loss of generality, we 
may assume an incoming wave that is left-incident on 
the diabatic state i = 1. Probability flow from Wi+ to 
the other CPWM components occurs only in the inter- 
action region, and then propagates outward towards the 
appropriate asymptotes. Adopting the usual normaliza- 
tion convention thus leads to the following boundary con- 
ditions: 

*i+(a;^-oo) = cxp +1 I y — <f>\ 

*(i>i)+(a!->-oo) = (13) 
-> oo) = 

In Eq. (fT3"]) above, is a time-dependent phase factor. 

As for the initial value, we take ^! a {x,t = 0) = 
for all a except a = (1,+), for which the basic WKB 
solution is used, as described in Ref. U Eq. (18). In the 
t — > oo limit, the resultant numerical solution \l/i±'s can 
be directly consulted to obtain transmission probabilities 
[associated with pi+(x — > oo)] and reflection probabilities 
[associated with pi-(x — > — oo)], as per Eq. ([T4")) . 

In most respects, the numerical algorithm employed 
for the time evolution is identical to that discussed in 
Sec. Ill Al One complicating factor is that the trajec- 
tories [determined from Vi(x)] are obviously completely 
different from one diabatic state, i, to the next. In a 
single-surface calculation, the + and — grid points move 
in opposite directions. Thus, the two grids must be in- 
commensurate for most times, although a simple scheme 
is used to ensure commensurate grids when t is a multiple 
of the shift time, t B hift- Starting at the left edge of the 
grid, xl, one propagates a single trajectory, x(t). The 
initial grid for both + and — is then taken to be the val- 
ues Xk = x{t = fcighift) for nonnegative integers, k such 
that Xk < xr, where xr is the right edge of the grid^ 

In the multisurface case, the trajectories at a given 
point x move with different speeds for different i values, 
thus rendering it impossible for grids to be commensurate 
across diabatic states, even at regular time intervals. A 
reasonable procedure, however, is to apply the above in 
component-wise fashion, i.e. by generating a family of 
/ trajectories, Xi(t). If the same i s hift value is used for 
all i, then at intervals t — fct s hift> the i+ and i— grids 
will be commensurate with each other, for all i (but still 
not across i values). The nonuniform grids generated in 
this fashion will be denser where u, is smaller — which 



is physically reasonable, and in any event is also true 
in the single-surface case. To evaluate the intersurface 
coupling contribution in Eq. ([7]), a polar intersurface grid 
interpolation scheme is used, exactly analogous to that 
used for (+<->—) coupling in the single-surface case. Note 
that all grids extend over the full coordinate range, i.e. 
roughly from xl to Xr. 

III. RESULTS 

In this section, we discuss the application of the nu- 
merical algorithms described in Sec. [II] to a variety of 
model systems. The mass m = 2000 a.u. was used in all 
cases. Also, fourth-order Runge-Kutta time integrators 
were employed, as was natural spline interpolation, un- 
less stated otherwise. For the symmetric, single-surface 
potentials, the phase-modified algorithm was used. A 
thorough and detailed convergence study was performed 
for each system, by varying each of the convergence pa- 
rameters in turn, and monitoring associated changes to 
the computed reflection and transmission probabilities, 
P z refl and p* rans . Unless stated otherwise, the precise con- 
vergence procedure followed is that described in Ref. 0, 
with one exception: since "oscillatory error" was always 
found to be negligible, reflection and transmission prob- 
abilities are defined here to be 

pr & = (^lMl)i*,-(^)i 2 

if-* = (v iR /v 1L )\* i+ (x R )\ 2 , (14) 
where v iL / R = lim Vi(x), 

x — >=foo 

rather than as described in Ref. |3- 

A. Single Surface Applications 

As fourth-order Runge-Kutta integrators have not 
been previously combined with the numerical innovations 
of Ref. 4J (Sec. Ill A"| . our first task is to evaluate the nu- 
merical efficacy of the present approach for single surface 
systems. In particular, we apply the method to a "test 
suite" of ID benchmark applications^ exhibiting a range 
of different attributes with regard to tunneling, barrier 
height and width, exoergicity, existence of reaction inter- 
mediates, and desired level of computational accuracy. 
In doing so, the robustness and stability of the method 
will be evaluated, as well as the numerical efficiency, as 
compared with the previous codes. 

Only a brief description of the individual test suite 
applications will be provided here; for additional details, 
the previous papers may be consulted. The five test suite 
potentials are as follows: 

• Eckart A: short, narrow Eckart barrier (height 
V « 0.0018 hartree) 

• Eckart B: tall, wide Eckart barrier (height Vq — 
0.011 hartree) 
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• Uphill ramp: tanh potential, barrierless and 
monotonic (Vr — Vl ~ 0.0018 hartree) 

• Barrier ramp: Eckart + tanh potential, asym- 
metric with barrier (Vr — Vl ~ 0.0018 hartree) 

• Double barrier: symmetric double barrier with 
reaction intermediate well (barrier height Vq rs 
0.0018 hartree) 

Results are presented in Table HI 

In the first investigation, we computed P refl and p trans 
for each of the five test suite potentials above, at an en- 
ergy E w V or E w (Vr- Vl) (Table[H Row 1) chosen to 
yield both substantial reflection and transmission. The 
results are presented in Columns II, IV, and VII-IX of 
Table H The computed P refl and P trans values (Table QJ 
Rows 8 and 9) were both converged to an absolute ac- 
curacy of around 10~ 4 (last digit uncertain). For the 
Eckart A system (Column II), all parameters were taken 
directly from Ref. d, except for the fixed (non-adaptive) 
Runge-Kutta time step size, A, which was reconverged. 
For all other applications, a full reconvergence of all pa- 
rameters was performed, as discussed above. For some 
systems, a direct comparison can be made between com- 
puted and exact (Rows 10 and 11) P rofl and P trans val- 
ues. For the Eckart systems, the exact results are known 
analytically^ For the uphill ramp system, "exact" val- 
ues were obtained from a much more accurate numerical 
calculation, converged to at least eight significant digits. 

In comparison with the corresponding calculations in 
Ref. 1, the converged parameter values are for the most 
part very similar, as expected. The glaring exception, of 
course, is the time step size, A, whose values are found 
here to be much larger than in Ref. 0, owing to the use 
of fourth-order time integrators. More significant, how- 
ever, is the fact that the present A values are substan- 
tially larger than those of Ref. H, which also employed 
non-adaptive fourth-order Runge-Kutta integration. For 
Eckart A for instance (Column II), the respective A val- 
ues (in atomic units) are as follows: Ref. |4|, 0.156; Ref. [3J, 
10.0; present work, 156. The latter value is indeed very 
large for the level of accuracy achieved, and can be at- 
tributed to the numerical improvements introduced in 
Ref. 13- Remarkably, an even larger A could in princi- 
ple have been used here, but the value given is already 
so large that A = t shift . Thus, over a single time step, 
each grid point moves to a neighboring site — resulting 
in an effectively "fixed" grid that obviates both spatial 
differentiation^ and interpolation. With regard to CPU 
effort, the Eckart A calculation performed here required 
4.5 ms on a 2.60 GHz Pentium CPU, as compared with 
1.7 s for Ref. 0. This reduction in CPU time is not as 
large as the 1000 x factor predicted by the increase in 
A, owing to the fact that each Runge-Kutta iteration re- 
quires additional computation, but is nevertheless very 
substantial; the present calculation is also orders of mag- 
nitude faster than the corresponding Ref. 3 calculation. 

In addition to the above investigation, we also exam- 
ined the present method's ability to provide extremely 



accurate results. This was done in two distinct ways. 
First, a very deep tunneling study of the Eckart B sys- 
tem was performed, as presented in Table U Columns V 
and VI. In effect, such a study evaluates absolute accu- 
racy rather than precision — e.g. Column VI, for which 
the transmission probability is around 10 -9 , but is com- 
puted to only three significant digits. The method clearly 
works well in this context — although the greater the ab- 
solute accuracy level, the less the improvement relative 
to Ref. 0- For the E = 0.1 Vq case, for instance, A in- 
creases by only a factor of 25 x, which is less than £ s hift) 
but only by a factor of two (Table Q] Row 6). Second, 
In order to evaluate precision rather than absolute ac- 
curacy, we also performed an extremely accurate calcu- 
lation of the E — Vo Eckart A system, for which both 
P refi and p trans were computed to eleven significant dig- 
its (Table [H Column III). As compared with Column II, 
the time step is reduced to an extent that is quantita- 
tively consistent with a fourth-order integration method. 
The corresponding increase in grid size is difficult to pre- 
dict quantitatively, but interestingly, is found to be such 
that the (£ s hift/A) ratio is again not far from one. Even 
for this extremely challenging application, the total CPU 
time required is only 36 seconds. 

Despite the substantial breadth of the test suite appli- 
cations considered, in all cases, the (t s hift/A) ratios are 
on the order of one (Table [H Row 6) — a situation that 
bears further scrutiny. Accordingly, we conducted an ad- 
ditional calculation of the Eckart A system using adaptive 
(Cash-Karp) fifth-order Runge-Kutta integration^ — the 
idea being that an adaptive method should automatically 
find the most appropriate time step sizes to use. In fact, 
using parameters similar to Column II, and specifying a 
target accuracy of 5 x 10~ 5 , we found the resultant A val- 
ues to be almost always equal to i s hift, except for the first 
few time steps, which had much smaller A values. The 
resultant computed transmission and reflection probabil- 
ity errors (relative to exact analytical values) conform to 
the target accuracy level. One might conclude from the 
above that the corresponding non-adaptive calculation 
using fixed A = £ s hift would yield similar results. In fact, 
the latter calculation results in a near identical p trans 
value, but a P rcfl with an error of 1.0 x 10~ 3 . Thus, the 
initially small time steps appear to play a vital role, and 
also argue cogently in favor of adaptive time integration 
schemes. 



B. Multisurface Applications 

We also performed numerical calculations for a vari- 
ety of benchmark multisurface applications for / = 2, as 
described below. 

• Pure coupling: asymptotically symmetric, with 
Vu(x) = and Gaussian coupling. 

• Tully model 1: crossed ramp potentials, Vu(x), 
with Gaussian coupling. 
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TABLE I: Parameters used for bipolar CPWM scattering calculations for single-surface test suite applications. All units are 
atomic units, except Row 12, the required CPU time in seconds on a 2.60 GHz Pentium CPU. Last digit of computed reflection 
and transmission probabilities (Rows 8 and 9) is uncertain. 



Quantity 






Sine 


le-surface test suite applications 








and 


Eckart A, E = V 




Eckart B, E = 




Uphill 


Barrier 


Double 


symbol 


low acc. 


high acc. 


V 


0.4 V 


0.1 V 


ramp 


ramp 


barrier 


energy. E 


0.001823 


0.001823 


0.011 


0.0044 


0.0011 


0.0023 


.0023 


0.0014 


grid size, N 


20 


800 


25 


61 


110 


19 


15 


20 


left edge, xl 


-2.0 


-4.0 


-2.6 


-3.5 


-3.5 


-1.5 


-1.5 


-2.2 


right edge, xr 


2.0 


4.0 


2.1 


4.0 


3.5 


2.2 


2.0 


2.2 


time step, A 


156 


2.47 


19.7 


14.9 


30.62 


107 


125 


196 


# steps, t shm /A 


1 


3 


3 


4 


2 


2 


2 


1 


max time, t max 


3899 


11867 


2893 


43978 


428621 


5792 


7000 


39143 


computed P rctl 


.28348 


.283.. .3842 


.45967 






.02391 


.45455 


.7958 


computed P trans 


.71665 


.716.. .6154 


.537 


1.560(-5) 


9.923(-10) 


.97595 


.54562 


.2053 


exact P rofl 


0.283358063869 


.459605 






.023901 






exact P trans 


0.716641936131 


.540395 


1.559(-5) 


9.920(-10) 


.976099 






CPU time (s) 


0.0045 


36 


0.035 


1.57 


15 


0.010 


0.008 


0.055 



• Tully model 2: double avoided crossing, with 
V22(x) a potential well. 

These systems exhibit a very broad range of attributes, 
especially with respect to energy scale. 

1. Pure coupling system 

The two-surface pure coupling system is defined via 
Vii(.t) = V 22 {x) = 

V 12 (x) = V 12 exp[-aOr - x )% (15) 

with the parameter choices Vq 2 = 150 cm -1 , xq — 0, and 
a = 1 a.u. This system represents an extreme case, in 
that there are no interaction potentials to induce di- 
rect (+<-»—) transitions. Instead, any reflection must 
arise indirectly, from a successive pair of intersurface 
transitions. As is clear from Eq. (|15p . the pure coupling 
system above is asympotically symmetric, implying that 
we can use uniform trajectories (which are also classical 
trajectories) and Eq. ©. 

The energy E = 100 cm" 1 = (2/3) V 12 is chosen to 
result in substantial probabilities for both scattered and 
unscattered wave components. Numerical calculations 
were performed using the algorithm of Sec. Ill C[ with 
adaptive Cash-Karp integration as discussed in Sec. lIII Al 
All reflection and transmission probabilities were con- 
verged to an accuracy of 10~ 5 or better, using the fol- 
lowing parameter values: N — 61; x^/r = =f3 a.u.; 
t max = 50000 a.u.; Runge Kutta tolerance level e = 10~ 6 . 

The four resultant converged component densities, 
Pi±{x) and P2±{x), are plotted in Fig. [1] as a function 
of x. All four curves are smooth and interference-free, 
as desired, and satisfy the appropriate boundary condi- 
tions of Eq. (fl~3|) . From Eq. (|14|) . the various reflection 
and transmission probabilities are found to be as fol- 
lows: P[ e& = 0.17886; P| cfl = 0.22382; Pf ans = 0.12194; 



T 




_l I I I 1 I L 

-3-2-10 1 2 3 

x (a.u.) 



FIG. 1: Component wave densities as a function of position, 
for the two-surface pure coupling system (Sec. IIII B ip . Cir- 
cles indicate trajectory grid points at the final time, for both 
diabatic states, 1 (filled circles) and 2 (open circles). 



ptrans = q.47537. In Fig. H the two transmitted wave 
component densities, p\+{x) and p 2 +{x), are summed 
together to form p + {x), and compared with the corre- 
sponding sum for the reflected components. As per the 
discussion at the end of Sec. Ill Bl the two summed curves 
are indeed found to be identical — apart from a constant 
shift, represented by the dot-dashed line. 



2. Tully Model 1 

The Tully models were introduced in the early 1990's 
by John Tully and coworkers, 32 to serve as benchmark 
numerical applications for investigating various processes 
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FIG. 2: Total transmitted (filled circles) and reflected (open 
circles) component wave densities (summed over all diabatic 
states) as a function of position, for the two-surface pure cou- 
pling system (Sec. IIII B ip . Circles indicate trajectory grid 
points at the final time, t — t m ax- Dot-dashed line indicates 
difference between the two curves — found to be a constant, 
equal to the total transmission probability. 



involving electronic transitions. 

The first Tully model is a simple avoided crossing sys- 
tem that consists of two diabatic ramp potentials — one 
an uphill ramp, and the other a downhill ramp — which 
cross each other symmetrically. For consistency and con- 
venience, we have replaced Tully's original piece-wise- 
cxponcntial ramps with an analogous tanh functional 
form. The modified Tully Model 1 potentials are there- 
fore as follows: 



V n (x) 

V 22 (x) 
V 12 (x) 



^ ~ 2 VL ^J tanl1 ~ 
^ Vl -Y^ j tanh [f3(x - x ) 
Vq 2 exp[-a(x - x ) 2 }, 



(16) 



with parameters Vr — —Vl — 0.01 hartree, x = 0, 
P = 1.2 a.u., Vq 2 — 0.005 hartree, and a — 1 a.u. 

The diagonal potentials, Vu{x), are no longer symmet- 
ric, necessitating the use of Eq. ([7]) rather than Eq. 
Due to the monotonicity of these potential curves, we 
take Vf s (x) = Vu(x), which results in standard classi- 
cal trajectories for the bipolar CPWM dynamics. Note 
that there are now potential couplings that can induce di- 
rect transitions between any pair of and VD^i com- 
ponents. However, the energy value chosen, E = 0.11 
hartree (« 24142 cm -1 ), is so much larger than the po- 
tential energy range as to ensure that there is negligible 
reflection. On the other hand, a substantial amount of 
intersurface coupling is obtained (Fig. [3]). 

Numerical calculations were performed as in 
Sec. IIII B 11 using the following parameter values: 




FIG. 3: Component transmitted wave densities as a function 
of position, and for a variety of times, t, for the two-surface 
Tully Model 1 system (Sec. IIIIB2|l . The upper family of 
curves represent pi+(x) at different times, whereas the lower 
family of curves represent p2+(x). The highest pi+ curve is 
the t — initial value (i.e. the WKB approximation), whereas 
the lowest p2+(x) = curve represents zero initial probability 
on diabatic state 2. Over time, probability is transferred from 
diabatic state 1 to 2 as indicated, via intersurface coupling. 



N = 50; x L/R = T3 a.u.; t max = 1000 a.u.; e = 10~ 6 . 
An "animation plot" of the resultant time-dependent 
transmitted wave densities, p\+{x,t) and p2+{x.i), is 
presented in Fig. [3l Each "contour line" is actually 
snapshot of a given component density at a particular 
time, t. For the uppermost curve is pi+(x) at the 
initial t — 0, whereas the lowermost pi+ curve represents 
the t — > oo limit. The opposite relation holds for the 
P2+ curves. From the figure, it is clear that the time 
evolution is smooth and well-behaved, and converges to 
the large-time limit quickly. The final curves represent 
essentially 100% transmission (the reflected wave den- 
sities are negligibly small, as expected), split roughly 
equally between the two curves as p* rans = 0.55016 and 
ptans = 0.44983. 

In addition to the E = 0.11 hartree case above, bipolar 
CPWM calculations were repeated and reconverged for 
a wide range of other energy values, to an accuracy of 
10~ 4 or better. Computed transmission probabilities for 
ptrans are presented in Fig.dJ and found to be in excellent 



agreement with the quantum results of Herman 



4!) 



3. Tully Model 2 

The second Tully model is a more complex, double 
avoided crossing system. The two diabatic curves consist 
of a symmetric well, and a constant energy curve which 
cuts across the well, giving rise to avoided crossings at 
the two intersections. These in turn give rise to quantum 



0.05 



0.1 0.15 
energy fa,U,j 



0.2 



0.25 



-3 



-2 

In(energy) 



-1 

;a.u.) 



E, for the two-surface Tully Model 1 system (Sec. IIII B 2]) . 
Filled triangles indicate present results, as obtained using the 
classical trajectory CPWM bipolar decomposition; solid curve 
indicates quantum results of Herman4£ 



interference effectsj 32 i 50 The potentials are 



Vu{x) 
V 22 (x) 



= 

= -V exp[-0(x - x Q ) 2 } + E 



y 12 exp[ 



-a(x — xqY 



(17) 



with parameters Vq = 0.10 hartree, xq = 0, f3 = 0.28 a.u., 
E = 0.05 hartree, V 12 = 0.015 hartree, and a = 0.06 
a.u. 

The diagonal potentials are symmetric, and therefore 
amenable to uniform trajectories. However, the energies 
involved are so large and "classical-like" that it is better 
to use standard classical trajectories, V° s (x) — Vu(x). 
Note that V 22 describes a potential well, thus ensuring 
that there will be no turning points. The large energies 
also ensure that reflection is negligible, but again, there 
is substantial intersurface coupling. Depending on the 
energy, the following range of parameters were used to 
achieve a convergence of 10~ 2 or better: N = 200-300; 
x L/R = =p8 a.u.; t max = 1000-2000 a.u.; e = 10~ 4 . 

Quantum interference associated with the two differ- 
ent surface pathways manifests as two distinct types of 
"Stiickelberg oscillations," 32 i 50 i5L52 h th associated with 
the Stiickelberg phase difference ) 52 ! 53 
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y/2m(E -Vu)- v / 2m(E - V 22 ) 



dx'. 
(18) 

Evaluating Eq. (fT8|) as a definite integral over the rele- 
vant range of x' yields oscillations in p* rans as a function 
of E. The computed transmission probabilities over a 
range of E values are presented in Fig. [5j these are in- 
deed oscillatory, and also in good agreement with Ref . H^. 

The second type of Stiickelberg oscillation is obtained 
by integrating Eq. {T5J to the indefinite limit, x, which 



FIG. 5: Computed transmission probabilities p^ rans vs loga- 
rithm of energy E, for the two-surface Tully Model 2 system 
(Sec. IIII B 3|) . Filled triangles indicate present results, as ob- 
tained using the classical trajectory CPWM bipolar decom- 
position; solid curve indicates quantum results of Herman. 4 - 
Note the Stiickelberg oscillations. 




FIG. 6: Component transmitted wave densities, pi+ (solid) 
and p2+ (dashed), as a function of position, for the two-surface 
Tully Model 2 system (Sec. IIII B 3]) . The oscillatory transfer 
of probability is a type of Stiickelberg phenomenon. 



can give rise to oscillations in pi(x) ~ pi+{x). Oscilla- 
tions are indeed observed in the converged transmitted 
wave densities for the E — exp(— 2) « 0.13533 hartree 
case considered, as presented in Fig. [6l The asymptotic 
oscillation wavelength as predicted by Eq. (|T5|) is found 
to be 1.31 a.u., in good agreement with the figure. 
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IV. CONCLUSION 

From a numerical perspective, the straight fourth- 
order Runge-Kutta bipolar CPWM algorithm, incorpo- 
rating the numerical refinements of Ref. 0, is found in 
Sec. MI Al to satisfy the demands of numerical stability, 
robustness, efficiency, and minimal user intervention. In 
particular, this algorithm is remarkably efficient at both 
high and low levels of absolute accuracy, and also for 
very high precision calculations such as Table [J Column 
III. We believe this algorithm to be the fastest avail- 
able for computing stationary scattering states of single- 
surface ID systems with predetermined boundary con- 
ditions, and advocate its use as a "black-box" method 
that can be applied to virtually any such application. 
That the algorithm also generalizes to multisurface ID 
applications in straightforward fashion (Sec. Ill Cj) . is also 
significant. We note that for systems such as Tully Model 
2, exhibiting substantial intersurface but little (+ «-> — ) 
interference effects, the bipolar CPWM in its present in- 
carnation necessarily undergoes a noticeable drop in ef- 
ficiency, as the Stuckelberg oscillations must manifest in 
the component field functions. However, Fig. [5] suggests 
that it may be possible to generalize the theory for ar- 
bitrary electronic representations (i.e. neither diabatic 
nor adiabatic) in such a way as to avoid interference os- 
cillations altogether in such cases. In any event, well- 
documented stand-alone fortran codes for the one- and 
two-surface algorithms considered here are available from 
the authors on request. 

Note that although little explicit connection with con- 
ventional "Bohmian mechanics" per se has been made 
in this work — i.e. of "polar" (amplitude-phase) decom- 
positions of the wavefunction components, and quan- 
tum potentials — implicitly, it is precisely this connection 
that enables the algorithms used here to work so well, 
as discussed in detail in the previous publications, no- 
tably Refs. 3 and 4. Of course, it is possible to recast 
Eqs. and © in polar form, yielding more explicitly 
"Bohm-like" results. Indeed, this is the usual convention 
in the prior semiclassical literature^ and we have also 
done this. However, this results in more complicated 
evolution equations that obscure aspects of the physics 
which we wish to emphasize in this work. Note that from 
a numerical perspective however, we use both polar and 
non-polar representations at different stages in the algo- 
rithm, transforming between them as described in Ref. 4 
[Eq. (9)]. 

From a theoretical perspective, the multisurface gen- 
eralization of the F approach with generalized classi- 
cal trajectories (Sec. Ill B[) presents some very intriguing 
aspects, with possible significance beyond the scope of 
the present paper. To begin with, the flux relations of 
Eqs. (fTQ|) through (|12| provide a natural, but non-trivial, 
multi-surface-like interpretation to the flow of probability 
among all 2f of the "components," ^i±. Even more com- 
pelling, however, is the form of the dynamical equations 



of Eq. ([9]), which above all serves to highlight the essen- 
tial sameness of the two types of transitions — intersurface 
and (+<-►—) — now treated together within a single uni- 
fied framework. That off-diagonal diabatic potentials in- 
duce transitions from one surface to another has always 
been known; now we find that in similar fashion, diagonal 
potentials also induce transitions, between transmitted 
and reflected components of the wavefunction associated 
with a given diabatic state. 

The above picture can lead one in a variety of in- 
teresting directions. Traditional TSH methods, for 
instance , 31 ! 32 utilize off-diagonal coupling 53 as a means 
of inducing trajectory hops from one surface to another. 
The present work suggests that a similar TSH scheme 
could be applied to effect transitions from reflected to 
transmitted wavepackets (even for single-surface calcu- 
lations), thereby naturally introducing interference and 
other effects that might otherwise be problematic in a 
traditional TSH context. Such a method would first re- 
quire that a non-stationary, wavepacket generalization of 
the present work be developed; indeed, this will serve as 
the focus of the next publication in the series, with the 
subsequent paper addressing multidimensional applica- 
tions. Note that both of these future papers will more 
directly emphasize the link with Bohmian mechanics — 
restoring the quantum potential, for example, which is 
required for wavepackets, but not essential for stationary 
state applications. As a bit of foreshadowing to motivate 
the present work, we comment that the multidimensional 
generalization still requires only two components per elec- 
tronic state (regardless of system dimensionality), and 
can accommodate standard Jacobi-type coordinate rep- 
resentations with arbitrary curvilinear reaction paths. 

Finally, we briefly mention two other possible areas 
for future development. First, the present theory is re- 
stricted to potentials that are asymptotically flat in both 
asympotes, x — > oo and x — ► — oc. It would be useful 
to generalize for potentials that are singular, or other- 
wise diverge in one asymptote or the other. Second, the 
intriguing result from Tabic U Row 6 that (i s hift/A) is 
always on the order of one suggests that simply setting 
this value to one would not adversely affect numerical ef- 
ficiency too severely. Indeed, such a modification would 
lead to some decided numerical advantages, such as the 
use of fixed grids which would avoid the need for intrasur- 
face (but not intersurface) interpolation. The resultant 
grid densities would be larger, however, and for this rea- 
son, such an approach might be substantially less efficient 
for high dimensionalities. 
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